%% Codes for the article "Beliefs and Financial Crises", by Arvind Krishnamurthy and Wenhao Li.
%% Set Path and Add Path (not needed if we call this file from "Main_Step1.m")
% mainpath = fileparts(mfilename('fullpath'));
% addpath(genpath([ mainpath, '/Codes_Model' ] ))
% cd(mainpath);

%% ------- Figure Appendix: Drift and Jump of Bayesian and Diagnostic Belief Process ---------------------------------------------------
S3 = load( ['solutions/baseline_behavioral.mat']);
sol_behavioral = post_processing(S3.model_sol);
par = sol_behavioral.par;
moments = simulation_results_output( S3, 1000);
% find the path of the lambda in the irrational case. 
dt = 1/1000;
lambda_grid = sol_behavioral.par.lambda_grid;
lambda_rational_sim1 = repelem( sol_behavioral.par.lambda_H, 10/dt)';  
for( iter = 1:(numel(lambda_rational_sim1)-1) )
    lambda = lambda_rational_sim1(iter);
    lambda_rational_sim1(iter+1) =  lambda +  sol_behavioral.par.mu_lambda_fun(lambda)*dt;
end
lambda_rational_sim2 = repelem( sol_behavioral.par.lambda_L*2, 10/dt)'; 
for( iter = 1:(numel(lambda_rational_sim2)-1) )
    lambda = lambda_rational_sim2(iter);
    lambda_rational_sim2(iter+1) =  lambda +  sol_behavioral.par.mu_lambda_fun(lambda)*dt;
end
lambda_rational_sim = [lambda_rational_sim1;  (lambda_rational_sim2) ];
ordering = [  1:length(lambda_rational_sim1),   length(lambda_rational_sim):-1:(length(lambda_rational_sim1)+1) ]';
% lambda_reference = [ lambda_rational_sim((end-1/dt):end),  lambda_rational_sim(1:(end-1/dt-1)) ];
cases_reference = ["avg", "low", "high"]; 
for(iter =1:numel(cases_reference))
    current_case = cases_reference(iter);
    if( current_case=="avg" )
        lambda_reference = mean(moments.lambda);
    elseif(current_case=="low" )
        lambda_reference = quantile(  moments.lambda, 0.1   );
    elseif(current_case=="high" )
        lambda_reference = quantile(  moments.lambda, 0.9   );
    end
    dambda_rational = sol_behavioral.par.mu_lambda_fun(lambda_rational_sim);
    lambda_behavioral_sim = irrational_lambda(lambda_rational_sim,  lambda_reference, 1, sol_behavioral.par );
    dlambda_behavioral = diff(lambda_behavioral_sim)/dt;   % drift of lambda in the behavioral belief case.
    dlambda_behavioral( length(lambda_rational_sim1) ) = 0;

    % Then also get the jumps of lambda, as a function of lambda_rational_sim.
    kappa_rational_lambda = sol_behavioral.par.kappa_lambda_fun(lambda_rational_sim);
    lambda_rational_sim_post_jump = lambda_rational_sim + kappa_rational_lambda;
    lambda_behavioral_sim_post_jump = irrational_lambda(lambda_rational_sim_post_jump,  lambda_reference, 1, sol_behavioral.par );
    kappa_irrational_lambda = lambda_behavioral_sim_post_jump - lambda_behavioral_sim;
   
    figure( 'Position', [0 0 360 300] ) ;  
    lambda_rational_vec = linspace( sol_behavioral.par.lambda_L,  sol_behavioral.par.lambda_H, 100 )';
    lambda_behavioral_vec = irrational_lambda(lambda_rational_vec,  lambda_reference, 1, sol_behavioral.par );
    plot(  lambda_rational_vec,lambda_rational_vec,  lambda_rational_vec, lambda_behavioral_vec ,    '--', 'LineWidth', 1.5 ); 
    hold on;   % xline(  min(lambda_rational_sim) );  text( min(lambda_rational_sim)+0.01, 0.4, 'minimum $\lambda^{Bayesian}$',  'Interpreter','latex' )
    xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{diagnostic}$', 'Interpreter','latex'); 
    saveTightFigure( strcat( figure_path ,'/model_illustration/behaviora_vs_rational_lambda_', current_case  ,'_reference_lambda.pdf' ) )

    figure( 'Position', [0 0 850 300] ) ;  
    subplot(1,2,1); plot( lambda_rational_sim(ordering), dambda_rational(ordering) ,  lambda_rational_sim(ordering(ordering<length(ordering))), dlambda_behavioral(ordering(ordering<length(ordering))),    '--', 'LineWidth', 1.5 )
    legend( {'rational belief', 'diagnostic belief'} ); xlabel('\lambda'); ylabel('$A(\lambda, \lambda_0)$', 'Interpreter','latex'); 
    subplot(1,2,2); plot( lambda_rational_sim(ordering), kappa_rational_lambda(ordering) ,  lambda_rational_sim(ordering), kappa_irrational_lambda(ordering),    '--', 'LineWidth', 1.5 )
    legend( {'rational belief', 'diagnostic belief'} ); xlabel('\lambda'); ylabel('$B(\lambda,  \lambda_0)$', 'Interpreter','latex'); 
    saveTightFigure( strcat(figure_path, '/model_illustration/lambda_process_illustration_with_', current_case  ,'_reference_lambda.pdf' ) )
end
close all;

%% -------- Figure Appendix: Diagnostic Belief as a Function of Rational Belief, with Differences in Reference Beliefs.
S3 = load( ['solutions/baseline_behavioral.mat']);
sol_behavioral = post_processing(S3.model_sol);
moments = simulation_results_output( S3, 1000);
% find the path of the lambda in the irrational case. 
dt = 1/1000;
lambda_grid = sol_behavioral.par.lambda_grid;
lambda_rational_sim1 = repelem( sol_behavioral.par.lambda_H, 8/dt);  % 10 years of simulation
t_vec = [1:numel(lambda_rational_sim1)]*dt;
for( iter = 1:(numel(lambda_rational_sim1)-1) )
    lambda = lambda_rational_sim1(iter);
    lambda_rational_sim1(iter+1) =  lambda +  sol_behavioral.par.mu_lambda_fun(lambda)*dt;
end
lambda_rational_sim2 = repelem( sol_behavioral.par.lambda_L*2, 8/dt); 
t_vec = [1:numel(lambda_rational_sim2)]*dt;
for( iter = 1:(numel(lambda_rational_sim2)-1) )
    lambda = lambda_rational_sim2(iter);
    lambda_rational_sim2(iter+1) =  lambda +  sol_behavioral.par.mu_lambda_fun(lambda)*dt;
end
lambda_rational_sim = [lambda_rational_sim1,  flip(lambda_rational_sim2) ];

% lambda_reference = [ lambda_rational_sim((end-1/dt):end),  lambda_rational_sim(1:(end-1/dt-1)) ];
cases_reference = ["avg", "low", "high"]; 
for(iter =1:numel(cases_reference))
    current_case = cases_reference(iter);
    if( current_case=="avg" )
        lambda_reference = mean(moments.lambda);
    elseif(current_case=="low" )
        lambda_reference = quantile(  moments.lambda, 0.1   );
    elseif(current_case=="high" )
        lambda_reference = quantile(  moments.lambda, 0.9   );
    end
    lambda_behavioral_sim = irrational_lambda(lambda_rational_sim,  lambda_reference, 1, sol_behavioral.par );
%     plot(t_vec,  lambda_rational_sim, t_vec, lambda_behavioral_sim)
    dlambda_behavioral = diff(lambda_behavioral_sim)/dt;   % drift of lambda in the behavioral belief case.
    dlambda_behavioral( length(lambda_rational_sim1) ) = 0;

    % Then also get the jumps of lambda, as a function of lambda_rational_sim.
    kappa_rational_lambda = sol_behavioral.par.kappa_lambda_fun(lambda_rational_sim);
    lambda_rational_sim_post_jump = lambda_rational_sim + kappa_rational_lambda;
    lambda_behavioral_sim_post_jump = irrational_lambda(lambda_rational_sim_post_jump,  lambda_reference, 1, sol_behavioral.par );
    kappa_irrational_lambda = lambda_behavioral_sim_post_jump - lambda_behavioral_sim;
    
    figure( 'Position', [0 0 360 300] ) ;  
    lambda_rational_vec = linspace( sol_behavioral.par.lambda_L,  sol_behavioral.par.lambda_H, 100 )';
    lambda_behavioral_vec = irrational_lambda(lambda_rational_vec,  lambda_reference, 1, sol_behavioral.par );
    plot(  lambda_rational_vec,lambda_rational_vec,  lambda_rational_vec, lambda_behavioral_vec ,    '--', 'LineWidth', 1.5 ); 
    hold on;   % xline(  min(lambda_rational_sim) );  text( min(lambda_rational_sim)+0.01, 0.4, 'minimum $\lambda^{Bayesian}$',  'Interpreter','latex' )
    xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{diagnostic}$', 'Interpreter','latex'); 
    saveTightFigure( strcat(figure_path, '/model_illustration/behaviora_vs_rational_lambda_', current_case  ,'_reference_lambda.pdf' ) )

end


%% -------------- Figure Appendix: Simulate the Model Under Different Variations of the Behavioral Belief --------------------------------------------------
moments = simulation_results_output( sol_behavioral, 1000 );
avg_lambda = mean(moments.lambda_rational);
simulation_key_moments_for_parameter_calibration;   % scripts that define the moments to be output, in the function "calibration_moments". 
lambda_vec = linspace( sol_behavioral.par.lambda_L,  sol_behavioral.par.lambda_H, 1000 )';
par = sol_behavioral.par;

% Plot the functional forms of each case. 
belief_fun1 = @(lambda) max( sol_behavioral.par.lambda_L + 0.001*lambda, min( par.lambda_H - 0.001*lambda,  lambda + 1.2 * sqrt( lambda - par.lambda_L ) .*  ( par.lambda_H - lambda )  )   );
belief_fun2 = @(lambda) max( sol_behavioral.par.lambda_L + 0.001*lambda, min( lambda,   2.46 * ( lambda - par.lambda_L ).^3  )   );
belief_fun3 = @(lambda) max( lambda, min( par.lambda_H - 0.001*lambda,  lambda + 2.3* sqrt( lambda - par.lambda_L ) .*  ( 2*avg_lambda - lambda )  )   );
belief_fun4 = @(lambda)  (lambda>avg_lambda).*min( par.lambda_H-0.001*par.lambda_H + 0.001*lambda ,lambda + 1* ( lambda -  avg_lambda ).^(1/3) .*  ( par.lambda_H - lambda )) + (lambda<=avg_lambda).* lambda   ;
belief_fun5 = @(lambda) min(lambda, 97* ( lambda - par.lambda_L ).^4)   ;
belief_fun6 = @(lambda) min( lambda, min( par.lambda_H - 0.001*lambda,  avg_lambda + 3.72 * ( lambda - avg_lambda ).^3  )   );

belies_list = { belief_fun1,  belief_fun2, belief_fun3, belief_fun4, belief_fun5, belief_fun6};
legend_list = {'pessimistic', 'optimistic', 'pessimistic', 'pessimistic', 'optimistic', 'optimistic'   };
for(iter = 1:length(belies_list))
    figure( 'Position', [0 0 360 300] ) ;  
    plot( lambda_vec, lambda_vec, lambda_vec, belies_list{iter}(lambda_vec),  '--' , 'LineWidth', 1); 
    xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{belief}$', 'Interpreter','latex'); 
    legend( {'rational', legend_list{iter}}, 'Location', 'northwest', 'Box', 'off' );
    saveTightFigure( strcat(figure_path, '/model_illustration/general_belief/case', num2str(iter) ,'.pdf') )
end
close all;  



%% Figure Appendix: Functional forms of Diﬀerent Beliefs Derived from Diagnostic Beliefs
lambda_reference = 0.2;
% Over-overtimistic case
S3 = load( ['solutions/baseline_behavioral.mat']);
model_sol = post_processing(S3.model_sol);
lambda_grid =  linspace(model_sol.par.lambda_L, model_sol.par.lambda_H,100 )';
lambda_rational = lambda_grid;
lambda_diag = model_sol.lambda_hat_fun(lambda_grid, lambda_reference);   
figure( 'Position', [0 0 300 300] ) ;  
plot( lambda_grid, lambda_rational,  lambda_grid, min(lambda_diag, lambda_rational) , '--', 'LineWidth', 1  ); xlim([0,0.7]);
lgd = legend('Bayesian belief', 'Overoptimistic belief' );  
xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{belief}$', 'Interpreter','latex'); 
% xlabel('Rational lambda (expected intensity of liquidity shocks)'); ylabel('Belief of liquidity shock intensity');
legend boxoff; lgd.FontSize = 10; lgd.Location='northwest';
set(gca,'FontSize',lgd.FontSize);
saveTightFigure( [figure_path, '/model_illustration/over_optimism_against_rational_lambda.pdf']  )

figure( 'Position', [0 0 300 300] ) ;  
plot( lambda_grid, lambda_rational,  lambda_grid,  lambda_diag, '--', 'LineWidth', 1  ); xlim([0,0.7]);
lgd = legend('Bayesian belief', 'Diagnostic belief' ); %  xlabel('Rational lambda (expected intensity of liquidity shocks)'); ylabel('Belief of liquidity shock intensity');
xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{belief}$', 'Interpreter','latex'); 
legend boxoff; lgd.FontSize = 10; lgd.Location='northwest';
set(gca,'FontSize',lgd.FontSize);
saveTightFigure( [figure_path,'/model_illustration/diagnostic_against_rational_lambda_fix_reference.pdf']  )

figure( 'Position', [0 0 300 300] ) ;  
plot( lambda_grid, lambda_rational,  lambda_grid,  max(lambda_diag,lambda_rational), '--', 'LineWidth', 1  ); xlim([0,0.7]);
lgd = legend('Bayesian belief', 'Overpessimistic belief' ); %  xlabel('Rational lambda (expected intensity of liquidity shocks)'); ylabel('Belief of liquidity shock intensity');
xlabel('$\lambda^{Bayesian}$', 'Interpreter','latex');  ylabel('$\lambda^{belief}$', 'Interpreter','latex'); 
legend boxoff; lgd.FontSize = 10; lgd.Location='northwest';
set(gca,'FontSize',lgd.FontSize);
saveTightFigure( [figure_path, '/model_illustration/over_pessimism_against_rational_lambda.pdf']  )



